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Transformation of the conventional radial Schrodinger equation defined on the interval r G [0, oo) 
into an equivalent form defined on the finite domain y(r) 6 [a, b] allows the s-wave scattering length 
a s to be exactly expressed in terms of a logarithmic derivative of the transformed wave function 
4>{y) at the outer boundary point y — b, which corresponds to r = oo. In particular, for an arbitrary 
interaction potential that dies off as fast as l/r n for n > 4, the modified wave function <f>(y) obtained 

by using the two-parameter mapping function r(y;f,/3) — f |l + -| tan(-7rj//2)J has no singularities, 

and 

2 1 d<t>(l)~ 



1 



7T/3 0(1) dy 

For a well bound potential with equilibrium distance r e , the optimal mapping parameters are f » r e 
and ft w § — 1- An outward integration procedure based on Johnson's log-derivative algorithm [B.R. 
Johnson, J. Comp. Phys., 13, 445 (1973)] combined with a Richardson extrapolation procedure is 
shown to readily yield high precision a s -values both for model Lennard- Jones (2n, n) potentials and 
for realistic published potentials for the Xe-e~, Cs2(a 3 Sj) and 3,4 He2(X 1 S^") systems. Use of this 
same transformed Schrodinger equation was previously shown [V.V. Meshkov et al. Phys. Rev. A, 
78, 052510 (2008)] to ensure the efficient calculation of all bound levels supported by a potential, 
including those lying extremely close to dissociation. 

PACS numbers: 31.15.-p; 34.50.-s; 37.10.De 



I. INTRODUCTION 

The so-called the s-wave scattering length a s is a key 
parameter for describing the interaction of particles at 
very low collision energies. In particular, the two-body 
collision problem is completely specified by the scatter- 
ing length in the low temperature limit where the elastic 
cross section becomes a e — 4™^ [l], 0. Many of the 
properties of a Bose-Einstein condensate Q also depend 
only on the scattering length. In particular, the chemical 
potential of a uniform Bose gas is simply proportional to 
the a s -value namely: fiBose — a s nA-n:h 2 /m, where n is 
the number density and m is the atomic mass Q. Thus, 
positive a s values correspond to overall repulsive interac- 
tions while negative values correspond to attractive ones. 

The scattering length is asymptotically related to the 
scattering phase shift r] s (k) by the expression [l|, H, @: 



= — lim 



tan i] s (k) 



(1) 



in which k is the relative wave vector of the colliding par- 
ticles. The scattering length can be also determined from 
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the wavefunction ip(r) of the radial Schrodinger equation 
d 2 ^(r) 



dr 2 
Q(r) 



2(i 



(2) 



n 2 



U(r) 



solved at zero energy E = with the inner boundary con- 
dition if>(0) — 0. It has been proved [l[ that for arbitrary 
potentials U(r) that obey the asymptotic condition: 



lim r n U(r) = for n > 3 



(3) 



the wavefunction tp(r) has the linear asymptotic form 



tp(r) ~ S (r - a s ) 



(4) 



in which the slope S is a constant. It can be seen from 
Eq.dU that the scattering length physically corresponds 
to the distance where the continuation of this asymptotic 
straight line crosses the r-axis. It is clear that the value 
of a s will be strongly depend on the interaction potential 
U(r) and on the reduced mass of the colliding particles 
fi. It is also well known that the scattering length a s is 
approximately related to the binding energy £^ max > 

nax) of the interaction 



of the highest bound level (v 
potential U(r) 0: 



H 2 



(5) 
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It is clear that a s —> +00 as E% — > 0, while a s will take 
on large negative values if the potential well is almost 
deep enough to support one more bound level. In other 
words, depending on the potential and the reduced mass, 
the scattering length can have any value on the interval 
(—00, +00). 

Accurate determination of scattering lengths deter- 
mined from photoassociation of ultracold colliding atoms 
8] offers the possibility of constructing reliable empiri- 
cal interatomic potentials all the way to the dissociation 
limit 0, ETJ]. However, rapid and robust methods for 
performing scattering length calculations are required to 
make such an inversion procedure feasible. Furthermore, 
a high degree of accuracy in the calculated a s -values is re- 
quired in order to provide the accurate derivatives of the 
scattering length with respect to the parameters defining 
the potential that are required for performing such fits 

n 

There are a number of existing schemes for calculating 
scattering lengths. However, calculation of the scattering 
phase shift r/ s (k) as a function of energy 0, Q in order 
to apply the low-energy extrapolation of Eq. (fTJ) leads to 
large uncertainties, especially for large a s values, since 
rj s (k) becomes a very steep function in the vicinity of 
k = 0. Moreover, the direct numerical integration of the 
radial equation ^ cannot be performed easily because 
the asymptotic linear form of Eq. (U|) is only reached at 
very large distances r, typically thousands of A. More ro- 
bust methods fllTfl3j are based on integration of Eq. @ 
to some finite matching distance r m , and then match- 
ing the resulting numerical wave function ip(r m ) with an 
asymptotic counterpart for the long-range region which 
is known analytically for particular kinds of potentials. 
Alternatively, the influence of a long-range interaction on 
the truncated apparent a s -value at some distance r = r m 
can be effectively corrected in the framework of a secular 
perturbation theory expansion [l4[. These asymptotic 
methods often yield fairly accurate results, but not for 
large a s values. Moreover, they all require additional 
computational effort to ensure the convergence of results 
with respect to an optimal r m value whose position is 
not well-defined a priori. These asymptotic methods 
also cannot readily be used for potentials for which the 
asymptotic solution is not available in closed form. We 
also note that a very elegant analytical formula for a s 
has been obtained using the semiclassical approximation 
for the wavefunction ij>(r) [15| . However, while that for- 
mula can be useful for estimation purposes, its accuracy 
is limited by its dependence on the WKB approximation 
[JEH, so ^ does n °t sumcc f° r many applications. 

The key to the present method is the analytically exact 
transformation fl7l EH of the initial Schrodinger equa- 
tion ([2]) defined on infinite domain r € [0, 00) into a 
modified radial equation defined by a reduced variable 
y(r) on a finite interval y(r) 6 [a, b] pjlHTj]. The result- 
ing transformed equation can be solved very efficiently 
by standard numerical methods [HI, Moreover, we 
show that a special choice of the mapping function r{y) 



leads to a transformed wave function <f>(b) that has no 
singularities on the interval [a, b]. In this case the scat- 
tering length can be expressed analytically in terms of the 
logarithmic derivative of the solution at the (finite) end 
point b and the associated mapping parameters. One ef- 
ficient way to obtain the required logarithmic derivatives 
of the wave function at y = b is to numerically integrate 
the resulting modified Ricatti equation using Johnson's 
log-derivative method [H, [24| and apply a Richardson 
extrapolation to the values obtained at the limit [1H Hil . 
A FORTRAN code applying this approach that has been 
applied both to model Lennard- Jones (2n, n) potentials 
and to realistic potentials for Cs2(a 3 E+) [12, El, EH, 
Xe-e" [13] and 3 ' 4 He 2 (X *£+) [H] taken from the litera- 
ture has been provided as Supplementary Documentation 
[29j | . An alternate version of the present approach based 
on the conventional Numerov propagator [21| has been 
implemented in a 'finite domain' version of the general 
purpose bound-state/Franck-Condon code LEVEL [30j . 

II. ADAPTIVE MAPPING PROCEDURE FOR 
SCATTERING LENGTH CALCULATIONS 

A. Reduced variable transformation of the radial 
equation 

We begin by introducing a mapping function y = y(r) 
that is a smooth, monotonically increasing function of 
the radial coordinate r and maps [0, 00) onto the finite 
domain [a, b\. The well-known substitution [13, EH 



^(r(y)) = M <t>{y) , g(y) = ^ > o (6) 

then transforms the conventional radial Schrodinger 
equation ^ into the equivalent form 

^ = - Q(v) m (7) 

in which 

Q(y) ee g 2 (y) Q(r(y)) + F(y) (8) 
whose additive term is defined as 

Hereafter, a prime (') symbol denotes differentiation with 
respect to our new radial variable y. 
It should be noted that: 

1. the modified wave function 4>(y) is now defined on 
the finite domain [a, b]; 

2. the modified radial equation ([7]) is completely 
equivalent to the initial one ©, as long as r(y) £ 
(C 3 [a, b] ) is a monotonically increasing function; 
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3. a knowledge of the inverse analytical function r(y) 
is all that is required to accomplish the exact trans- 
formations of Eqs. ©-((9]); 

4. for a fixed mesh of equally spaced y points on the 
interval y € [a,b], the function p(r) — dy/dr = 
l/5(y) defines the density distribution of mesh 
points of the initial coordinate r. 



It is clear that conventional finite-difference 121, 23, 2 



|26| ] and pseudospectral [22j methods can be used straight- 
forwardly for integrating the transformed radial equation 
out to the point y — 6, as long as the modified wave 
function <p(y) is not singular anywhere on the interval 
[a, b], and especially not at its end points. Furthermore, 
it will be proved in the next section that the scattering 
length a s is an explicit function of logarithmic derivative 
of the wave function <p(y) 



ay) 



<Ky) 



(10) 



at the boundary y — b. It is well-known that the trans- 
formation (|10p converts the radial equation ([7} into the 
modified Ricatti equation 



e(y) + Q(y) + e(y) = o 



(ii) 



which can be numerically integrated, for instance, using 
Johnson's efficient log-derivative method [13, HiJ • More 
details regarding this method are given in the Appendix. 



B. Developing a Formula for the Scattering Length 

Since Eq. (Q| shows that the scattering length can be 
expressed formally in terms of the asymptotic behavior 
of the ordinary radial wavefunction ip(r), we can write 



dijj(r)/dr 



(12) 



It therefore seems desirable to investigate the asymptotic 
behavior of the modified function 4>{y) near the outer 
boundary point y = b corresponding to r — > 00, where 



0(5) 



S 



r(y) 



y/dr(y)/dy 



y^b 



(13) 



Since we are mapping an infinite interval onto a finite 
one, it is clear that the mapping function r(y) must have 
a singularity at the upper end of the interval. To proceed, 
we assume that this singularity has the form 



1 



(b - vY 



y 



(14) 



where necessarily 7 > 0, because the r{y) function must 
go to infinity as y approaches b in order to convert the fi- 
nite domain y £ [a, b] into the infinite one r £ [rmin, +00]. 



Inserting Eq. ([14} into Eq. (U3]| yields 



4>{y) ~ S [b-y) 1 ^ 



y^b 



(15) 



and differentiating that result with respect to y yields 
!-7 



4>'{y) 



s 



(b-y) 1 



(16) 



From Eqs. fTS]) and ([TB]) it is immediately clear that for 
both the modified wave function 4>(y) and its derivative 
4>' (y) to be non-singular at the point y = b (or r = +00), 
necessarily 7 = 1. 

Avoidance of this singularity at y = b is clearly desir- 
able if we are to integrate the modified radial equation 
(0 accurately on the whole interval [a, b], so it is nec- 
essary to choose a mapping function r(y) that has only 
a simple pole at the end point y = b. Let us assume 
that r{y) can be expressed as a Laurent series expansion 
about this point: 



r(v) 



C-l 

b-y 



CO 



X>(6-y)* 



(17) 



where c_i > 0. Inserting Eq. (TTT)) into Eq. (fTB")) and dif- 
ferentiating the resulting expression for <p(y) with respect 
to y yields, at the point y = b, 



which shows that 



0'(6) 



S(a s - c ) 



C() 



,(18) 



(19) 



is the desired relationship between the scattering length 
a s and the log- derivative function at the outer end of the 
interval, £(&). 

To facilitate accurate calculation of the the required 
£(&) value, it is desirable that the function Q(y) of Eq. dSJ 
be nonsingular at the end point y — b. It is easy to verify 
that the F(y) contribution is nonsingular at y = 6, since 



F(b) 



3ci 

C-l 



(20) 



Moreover, from Eq. (IT71) it is clear that in the limit y —> b 

(or r — ¥ 00), 



g(y) 



dr{y) 
dy 



C-l 



(b-y) 2 



oc r 



(21) 



This means that for any potential which dies off more 
slowly than 1/r 4 , the product g 2 U(r) that comprises the 
main part of the function Q(y) (see Eq. ©) diverges in 
the limit y — > b (r — > +00). In particular, if U(r — > 
00) ~ — C n /r n , then as y — > b , 



9 2 U{y) * £ 



2/i CJb-yY 



h 2 (c_i)' 



(22) 
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and hence 



Q(b) = 



F(b) 

H 2 ( C _!) 2 
DO 



n > 4 
F(b) n = A 



This means that, to calculate the modified wavefunction 
<p(y) and/or its derivative (f>'(y) at the end of the range 
(y = b) for an n ~ 4 case requires explicit knowledge 
of the leading long-range induction coefficient C4, while 
for the more common n — 5 or 6 cases it requires only 
a knowledge of the mapping function r(y) that defined 
F(b). In contrast, for n < 4 no satisfactory determi- 
nation of </>(&) and/or 4>'(b) can be achieved using the 
present approach. This latter result is consonant with 
the fact that the scattering length is not defined for po- 
tentials that die off more slowly than 1 /r 4 [l| . 



C. Introduction of Mapping Functions 



To make practical use of our relationship (119)) between 
the scattering length and the log-derivative function, it 
is necessary to introduce an analytical mapping function 
with the Laurent expansion form of Eq. (|17[) . The sim- 
plest way to do this is to define the mapping function as 
the first two terms of the Laurent series: 



r(y) = 



y 



C{) 



(23) 



since only the c_i and cq coefficients are required in 
Eq. (fT9l) . In particular, setting 6=1 and c_i/2 = — c = 
r, in ([23)) yields a version of (omitting a factor of 2) of 
the well-known Ogilvic-Tipping (OT) potential energy ex- 
pansion variable [1 91 ] : 



_ T — V 

Voiir; f) = — — ; y Q1 £ [-1, 1] 



(24) 



For this case Eq. (fl9|) shows that the scattering length is 
explicitly defined as 

o, = r[2£(l)-l] (25) 

while the reciprocal mapping function is 

'l + y s 



roi(y) 



and 



i-y 



2r _ (r + f) 2 



2f 



(26) 



(27) 



Moreover, for this mapping function the additive term 
F OT (y)=0. 

Use of the simple one-parameter mapping function of 
Eq. (|24p does provide a reliable method for calculating a s . 
However, our recent experience with applying this type 



of transformation to bound-state problems [l8| suggests 
that better efficiency may be attained using more flexible 
two-parameter functions. One such function which was 
successfully used for bound-state calculations is the two- 
parameters (r, a > 0) Surkus variable [20|: 



ys(r;r,a) = 



j.UL _l_ jtUL 

for which the reciprocal mapping function is 



ysr€ [-1,1] (28) 



rs(y) = r 



1 



y 



i-y 



l/a 



(29) 



However, except for the special case a = 1 in which 
yg(r;f, a) reduces to ym{r]f), this mapping function is 
not appropriate for scattering-length calculations, since 
when a / 1 the mapping function rs(y) does not take 
on the required Laurent series expansion form (|17[) as 
y — > 1, so the associated log-derivative wavefunction 
becomes singular there. 

In contrast, the two-parameters tangential function 

m 



rt s {y) 



1 



1 /ny 



(30) 



defined on the interval y £ [(2/7r)tan 1 (— /3),1], whose 
inverse is the inverse tangent function 



tan" 1 [(3 (r/r- 1)] 



(31) 



completely conforms to the required Laurent expansion 
form of Eg. (fT7|) as y 1, since 



r t g (y) 

For this case 



2f 



- y) 
"2^(1 



7T/3 



7TT 

t (1 — y) 

6/3 1 J> 



1 



(32) 



(33) 



and it is easy to show that the additive contribution to 
Eq. (jHJ is a constant, 



7T - 

T 



and that 



7TT 

3tg(y) = ^ cos" 



(wy/2) 



(34) 



(35) 



nr 

2ft 



1 + ft 2 (r/f- iy 



D. Determination of Optimal Mapping Parameters 

The optimal values of the parameters defining the map- 
ping function of Eq. (|30p may be expected to depend both 
on the nature of the interaction potential U(r) and on 
the particular numerical method used for integration of 
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FIG. 1. Comparison of the original ip(r) and modified </>(y) = 
ip/^/gig unbound wave functions at zero energy for the a 3 St 
state of Cs2 [H| , as calculated using the Numerov method [2ll . 
|29| ]. The scaling function </t g (7*) was calculated from Eq. (I35p 
using the mapping parameters f = r e = 12.0 [ait] and j3 = 2. 



the modified radial equation ([7]). When u sing Johnson's 
log-derivative method of integration [23l . l24j . the opti- 
mal mapping parameters r and /3 could be determined 
by minimizing the truncation error A£ estimated using a 
Richardson extrapolation (RE) (|A.7[) . or by minimizing 
the difference between (y = b) and £h 2 (y — b) values 
corresponding to two different integration step sizes hi 
and h%, 



min \€hi ~ 6i 2 



(36) 



Note that it is tacitly assumed that the optimal param- 
eter values do not depend on the integration step size. 
This has been confirmed by the illustrative calculations 
presented in Section Mil It will also be shown there that 
the minimum of the functional (|3l)]) is a rather smooth 
function of both mapping parameters, f and j3. 

It seems reasonable to assume that the optimal map- 
ping function would correspond to the situation in which 

Q(y) = g 2 (y) Q(r(y)) + F(y) « const (37) 

since it is well-known that minimal truncation errors arise 
in particle- in-a-box problems where the potential U(r) 
does not depend on r. Furthermore, within the classically 
allowed region where Q>0, the additive term F(y) can 
be approximately neglected [l|, [l6[ . In this case Eq. (|3"T)l 
reduces to g 2 Q const, and inverting this condition 
(recalling that g(y) = dr/dy) yields 



y op t{r > r ) 



(38) 
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FIG. 2. (Color online) Mapping density functions p(r) — 
dy/dr = l/g(r) for the a 3 E+ state of CS2 [la] calculated: (i) 
within the frame work o f the conventional WKB approxima- 
tion pwkb = \l —U (r) p9p . (ii) using the one-parameter 
ym(r;f) mapping function (|24p with f = r e = 12.0 [aw], 
and (Hi) using the two-parameters ytg(r; f, P) mapping func- 
tion (|31|) with p = r e and /3 = 2. i/»(r) is the associated 
zero-energy wavefunction. 



Here ro is the left-hand (inner) classical turning point, 
which is the root of the equation U(r ) = 0. 
Differentiating Eq. (|38|) with respect to r yields 



Popt 



dy a 



P t 



dr 



VQ{r 



%V(r) (39) 



For optimal mapping, therefore, the density function 
Popt(^) should be close to Q(r) in the Q(r) ^0 region 
where the original wavefunction ij)(r) has its maximum 
oscillation frequency. From Eq. (|3^|) it follows that the 
optimal mapping does not depend on the reduced mass 
jj,, because it appear in Q(r) simply as a multiplicative 
factor. It is also easy to see that the mapping function 
(|38p transforms the modified wavefunction <fi(y) into the 
familiar particle-in-a-box form 



<£o P t(y) ~ ^sin(fcy) + Bcos(ky) 



(40) 



which correlates with the conventional WKB approxi- 
mation [l|, [lj| in the classically allowed region where 
Q(r) > 0. As is shown by Fig. 1, in contrast to the 
original wavefunction ij)(r), the modified wavefunction 
4>(y) has loops of almost constant amplitude and spacing 
over the potential well, behavior which is qualitatively 
very similar to that implied by Eq. (|4"0")l . It is therefore 
expected that the WKB mapping of Eq. (38) should be 
close to "optimal" for all numerical methods (such as the 
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finite-difference and collocation methods) that are based 
on an equidistant grid in the classical region of motion. It 
should be stressed, however, that the mapping function 
defined by Eq. (|38p cannot itself be applied for scatter- 
ing length calculations because it does not satisfy the 
required asymptotic behavior of Eq. (fl~7|). and because 
it becomes imaginary in the classically forbidden region 
where r < tq. 

It is easy seen from Eqs. (|27p and (1551) that the den- 
sity function pm(r) corresponding to the one-parameter 
mapping function of Eq. f|24[) is proportional to l/[r + r] 2 , 
while that for the two-parameters tangent function of 
Eq. (|31[) has a Lorentzian form with a maximum at f. 
The plots presented in Fig. [51 show that there is good 
agreement between the pwkb(?~) and ptg(^) density func- 
tions over much of the domain of the wavef unction i/>(r). 
This leads to the conclusion that the optimal value of 
parameter f opt should be close to the equilibrium inter- 
nuclear distance r e of the potential U(r). Hence, for the 
two- parameters mapping function (|3T|) . it seems reason- 
ably to fix f op t ~ r e and to vary only the single parameter 
p. Note that the pronounced differences among the den- 
sity functions for these three cases at small distances is 
not very important, because the associated wave function 
ip(r) dies off exponentially in this region. 

In the following, therefore, parameter f of the mapping 
functions ym{r] f) of Eq. (|24|) and yt g (r; f, (3) of Eq. ([3"T]) is 
fixed as f = r e , while the optimum value of /3 of Eq. (|3"Tj) 
is determined from a one-dimensional minimization of the 
functional of Eq. (f36| . 



III. 



IMPLEMENTATION, TESTING, AND 
DISCUSSION 



A. Model-Potential Applications and Tests 

A computer program based on the adaptive mapping 
procedure described above has been written and tested 
29] , both for a variety of model potentials, and on poten- 
tials for real systems taken from the literature. All results 
were obtained on 32 bit processors, mainly using double- 
precision arithmetic, but with quadruple-precision calcu- 
lations being used sometimes in order to delineate the im- 
pact of accumulated numerical round-off error. This sec- 
tion presents the results of those illustrative applications, 
together with some general discussion of the method. 

The first sample calculations presented here are for the 
model Lennard-Jones(2n, n) potentials 



TABLE I. The s-wave scattering wavelengths a s (in A) calcu- 
lated for the LJ(2n, n) potentials defined by Eq. (f4"Tj) [29| . All 
models have equilibrium distance r e = 1 [A] and the reduced 
mass is set as /i = 16.85762920 [au] so that in "spectroscopists 
units", the scaling factor h ' 'fan = 1 [cm -1 A ]. ® e is the well 
depth in cm -1 , while Umax is the vibrational quantum number 
of the last bound level supported by the potential. 
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energy approximation of Eq.©. 
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FIG. 3. (Color online) Convergence tests for a s -values cal- 
culated in double precision arithmetic for the model 15-level 
LJ(2n, n) potentials of Table U using the one-parameter map- 
ping function ym(r\f) (|24[) with f = r e = 1 [A]. Calculations 
performed using Johnson's log-derivative method alone are 
denoted 'JLD', while the results labeled 'JLD-RE' were ob- 
tained by also applying a (A,A/2)-Richardson extrapolation 
procedure to those results. 



^Lj(r) 



1), 



(41) 



defined by the potential function parameters and system 
reduced mass listed in Table HI The first three of these 
LJ(2n,n) potentials, each supporting 15 bound vibra- 
tional levels, are the same models systems considered in 
our recent application of this adaptive mapping approach 



to bound-state problems [Tj|. The last two LJ(12, 6) po- 
tentials have much deeper wells and support 100 vibra- 
tional levels. They were introduced here to highlight the 
efficiency of the present method, since the computational 
effort of scattering-length calculations increases dramat- 
ically as D m ax and the absolute magnitude of a s -values 
increase. 
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FIG. 4. (Color online) Convergence tests for a s -values cal- 
culated by the JLD-RE(A r ,A r /2) method in double precision 
arithmetic (solid curves) for the five model LJ(2n,n) poten- 
tials of Table U using the two-parameters ytg(r;r = r e ,P) 
mapping function of Eq. ()31[l with the optimal parameters 
/3 op t listed there. The dashed lines present results obtained 
using quadruple-precision arithmetic. 



Converged a s -values for five model LJ(2n,n) poten- 
tials are presented in the last column of Table HI As 
discussed above, the range-mapping parameter was fixed 
at f = r e = 1[A] for both yo-r(r) and ytg(i~) mapping func- 
tions. The optimal /3 op t values for the latter case were de- 
termined by minimizing the functional of Eq. (|36p , yield- 
ing the results shown in the column four. It is interesting 
to see that these empirically determined (3 opt values have 
essentially the same dependence on the inverse-power n 
governing the long-range behavior of the potential as was 
the case for the analogous parameter a of the Surkus- 
variable mapping ([28]) used in the bound-state study of 
Ref.Q. 

The log-log plots in Figs,[3] and |4] show how the relative 
errors in the calculated a, values 



S = 



-.calc 



(42) 



depend on the number of mesh points N used in the 
radial integration. The reference values a!! xact were ob- 
tained from calculations using quadruple precision arith- 
metic by increasing N until the full desired level of con- 
vergence was achieved. The numerical 'noise' on these 
plots for S < 10~ 12 indicates the precision limits achieved 
using ordinary double-precision arithmetic, while the 



dashed lines show how the convergence trend continues 
when quadruple-precision arithmetic is used. 

The three upper curves in Fig. [3] (labeled JLD) display 
results obtained using Johnson's log-derivative method, 
while the three lower curves illustrate the greatly im- 
proved convergence achieved when that approach is cou- 
pled to the (JV, JV/2) Richardson extrapolation proce- 
dure described in the Appendix. These results clearly 
confirm the prediction of Eq. (IA.5|) . that Johnson's log- 
derivative method (JLD), and that method combined 
with a Richardson extrapolation to zero step ( JLD-RE) , 
demonstrate A^ -4 and TV -6 convergence rates, respec- 
tively. We see that for these 15-level LJ potentials, the 
JLD-RE(Af, N/2) method allows us to attain 11-12 sig- 
nificant digits in calculated a s values when using only 
N « 10 4 radial mesh points. 

The results displayed in Fig. [4] show that as expected, 
the computational effort increases significantly for the 
100-level LJ(2n, n) models, especially for the very last 
case considered in Table 1 , for which the scattering length 
is extremely large: a s ss 10 4 [A]. However, that even for 
those cases full double-precision converge is achieved with 
only N « I0 5 . Note that for a given number of grid 
points, use of the simple one-parameter mapping func- 
tion yar(r) yields relative errors that are only ~ 10 times 
larger than the analogous results yielded by the optimum 
two-parameter functions yt g (r). The non-linear 'cusp'- 
like behavior near A~ rj 10 3 on the JDL-RE plots in 
Figs. |31 and 0] appears to be due to accidental cancela- 
tion of higher-order terms in the RE series expansion for 
these cases. It does not appear in analogous results based 
on use of a Numerov propagator. 

The results presented in Fig. [5] clearly show that the 
accuracies of a s values calculated with any given num- 
ber of mesh points A~ vary smoothly with the value of 
mapping parameter parameter /?. They also show that 
the optimal parameter value /3 opt is essentially indepen- 
dent of the number of integration points used. Analogous 
results are obtained when modeling LJ(2n, n) potentials 
with n = 4 and 5. 

The sensitivity of the value of /3 pt to the depth of the 
potential energy well has also been investigated. Figure[S] 
shows that the values of /3 op t for model LJ(2rt, n) poten- 
tials only depend significantly on the dissociation energy 
S) e for very shallow potentials which support only a few 
bound states. For such species (e.g., see the He2 ex- 
ample considered below), accurate scattering lengths can 
readily be calculated using a relatively small number of 
grid points, so determining precise values of /3 pt is im- 
material. Overall, Fig. [5] shows that the optimal value 
of the /3 parameter in the mapping function of Eq. (|3 1[) 
depends mainly on the inverse-power n governing the lim- 
iting long-range behavior of the interaction potential. As 
indicated by Fig. [5] and Table HI these optimal values are 
approximately defined by the simple relationship 

--1 , 
2 



-"opt 



(43) 



at least for these LJ(2n, n)-like potentials. We believe 
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FIG. 5. (Color online) Relative errors (|42|l in the a s -values 
obtained for the 15-level LJ(12, 6) potential (|41|l as functions 
of the mapping parameter /3, where A*' is the number of grid 
points used in the JLD method. 
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FIG. 6. (Color online) Optimal mapping parameter /3 op t as 
a function of the dissociation energy 55 e of 15-level LJ(2n, n) 
potentials (|41[) for n = 4, 5 and 6. 



that this empirical rule will be valid when applying 
the two-parameter mapping function of Eq. (|30[) to any 
deeply bound interatomic potential having a long-range 
tail which dies off as 1/r™ for n > 4. 

It is interesting to note that the two-parameters map- 
ping function yt g (r;f,f3) of Eg. (|31l) is also well suited 
for solving the Schrodinger equation for the bound-state 



FIG. 7. (Color online) Optimal mapping parameters /? op t de- 
termined for the bound vibrational levels v G [0, 14] of the 
model 15-level LJ(2n,n) n=4,5,6 potentials (|41[) . The dashed 
horizontal lines corresponds to the /3 op t values implied by 
Eq.flSJ). 



problem discussed in Ref. In particular, Fig. [7] shows 
how the values of (3 op t evolve if this mapping parameter 
is optimized independently for each level of our three 
15-level LJ(2n, n) potentials. Except for the very last 
level, the values of /3 pt smoothly approach the limiting 
value implied by Eq. ([4"3"f as the vibrational levels ap- 
proach dissociation. The substantial deviation from this 
limiting behavior for the lower vibrational levels is not 
a matter of concern, since those levels can readily be lo- 
cated accurately using a very wide range of f3 values, and 
since their radial amplitude is much smaller than that for 
the highest levels, the overall computational efficiency for 
such levels is not very strongly dependent on the value 
of f3. The abrupt drop-off in the values of /3 op t for the 
very last level also mimics the behavior found in Ref. [l8[ 
using very different mapping functions (|29[) based on the 
Surkus variable [2(| defined by Eq. (l2"51) . This abrupt de- 
crease of /3 op t for the last very weakly bound vibrational 
level v = i>max = 14 is related to the form of the cor- 
responding wavefunction ipv m ^(r) which has very broad 
last loop centered at a very large internuclear distance. 
In particular, for the case of a last level which lies ex- 
tremely close to dissociation, it appears that a small lim- 
iting value of /3 op t < 0.5 is needed to provide a sufficiently 
broad distribution of grid points to properly characterize 
the outermost loop of the wavefunction. 
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B. Applications to 'Real' Systems 

This section describes application of our new method 
of performing scattering length calculations to three 'real' 
physical problems. The first of these is the elastic scat- 
tering of a free electron from a neutral Xe atom. The 
function used for the Xe-e - interaction potential is the 
'HFD'-type potential reported by Szmytkowski [l3| (with 
energy and distance in atomic units): 

EWM = Ae^ r2 - F n (r)C n /r n , (44) 

n=4,6 

in which A = 306.0 and 7 = 1.0, while the n-dependent 
damping function F n has the form 



F n {r) 



1 



-(r/r c y 



(45) 



with r c = 1.89 being a cut-off radius, such that F n (r) = 
for r < r c . Following Czuchaj et al. its disper- 

sion coefficients C4 = «i/2 and Cq — (a?, — 6/?i)/2 in 
Eq. (|44l correspond to the familiar charge/induced-dipolc 
and charge/induced-quadrupole interactions, where a\ = 
27.292 and «2 = 128.255 are the static dipole and 
quadrupole polarizabilities of the Xe atom, while f3\ = 
29.2 is the dynamical correction to the dipole polariz- 
ability. This shallow U Xe . e - (r) interaction potential sup- 
ports no bound levels and its scattering length is nega- 
tive, as was determined by Szmytkowski using an asymp- 
totic method [T3 |. 

The present a s -value for this system (in au), given in 
the first row of Table Hi) was converged to 14 significant 
digits using the JLD-RE(7V,iV/2) method based on N « 
10 3 grid points and a scaling factor of H 2 /2/x = 0.5. How- 
ever, only N 100 grid points are required to converge 
the present a s -value to 7 significant digits. The present 
estimate agrees with the value reported by Szmytkowski 
13] to within the 6 significant digit that he reports (see 
row 2 of Table HU). 

Our second real-world example is a modified 'Hartrcc- 
Fock dispersion ' (HFD) type potential for the a 3 S+ 
ground triplet state of the cesium dimer [12, [H, UH 
(again, with energies and distances in atomic units): 



U Cs2 (r) = Ar a . 



F(r) 



E 

n=6,8,10 



C n /r n 



(46) 



The first term in Eq. (|46j) . defined by constants A = 8 



10" 



5.53 and 7 = 1.072, represents the exchange 



repulsion energy, while the second is a sum of van der 
Waals dispersion terms (with coefficients Cq — 7.02 x 10 3 , 
C 8 = 1.1 x 10 6 , C10 = 1.7 x 10 s ) multiplied by the n- 
independent damping function 



F(r) = H(r- 



H(r c - r) e 



-(r c /r 



(47) 



where H(x) is the Heaviside step function: H(x) = 1(0), 
when x > (<)0, while r c — 23.165 is the cut-off ra- 
dius. This potential energy function supports up to 58 



TABLE II. Comparison of s-wave scattering lengths a s (in 
ait) obtained in the framework of JLD-RE(iV,iV/2) proce- 
dure [2S| using realistic published potentials for Xe-e - [27|, 
Cs 2 (a 3 Ej) pSGlGl, and 3 ' 4 He 2 (X 1 S+) [H interactions 
and two-parameter mapping function yt s (r;f,/3) of Eq. (|31[) . 
The optimal mapping parameters /3 op t were determined by 
minimization the functional (|36[) . while the f values (in au) 
were fixed at f op t = r e , where r e is the equilibrium distance 
of the relevant potential. 



species potential n v n 



Popt I opt 



Xe-e" Eq. dm 4 1 -4.95280521509712 1.7 3.2 

-4.9528r 

Cs 2 Eq. (gUl 6 57 68.215967213 2.0 12.0 

68.21596^. 
68.21823~ 
68.0^. 

4 He 2 Eq. gSJ 6 236.3688418603 2.5 5.67 

236.36884174 d 



3 He 2 



* -13.17845206095 3.6 5.67 
-13.178452062 d 



There are not bound levels. 



a Ref . 13 
b Ref. IS 
c Ref. 15 
d Ref. E 



bound levels, and can be considered a typical example of 
a many-level interatomic potential. 

The present et s -value, obtained using our JLD- 
RE(iV ,N/2) procedure with a scaling factor of H 2 /2/j, = 
1/2.422 x 10 5 , is listed in row 3 of Table HI It agrees 
with previous estimates obtained using WKB [l5| (row 
6), asymptotic [l3[ (row 5), and iterative [l2[ (row 4) 
methods to within 2, 4 and 7 significant digits, respec- 
tively. The present method required about N m 10 4 grid 
points to obtain 10 significant digits in the calculated 
a s -value. 

Our third practical application is to the ground X 1 E+ 
state of the helium dimer, for which the interaction po- 
tential is again represented by an 'HFD'-type function 
[28| written as (with energies in K and distances in A) 



U He2 {r) =2) 



Aexp(— jx) 



. (48) 



in which x = r/r m , D = 10.8, A = 544850.4, 7 = 
13.353384, C 6 = 1.3732412, C 8 = 0.4253785, C 10 = 
0.1781, and r rn — 2.9673. The damping function F(x) 
in Eq. (|48p is defined by Eq. (|4T)) with the parameter 
x c = r c /r rn = 1.28. This is an exotic example of very 
shallow interatomic potential, as it supports only a sin- 
gle bound level for the heavier isotopologue 4 He2 and no 
bound levels at all for the lighter isotopologue 3 He2 . For 
consistency with the calculations of Ref. ll| . the values 
of the scaling factors fi 2 /2/x used in the present calcula- 
tion were 16.085775 and 12.120904 for 3 He 2 and 4 He 2 , 
respectively. 
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FIG. 8. (Color online) Comparison of convergence behav- 
ior for scattering lengths implied by the 3 Hc2(A' poten- 



tial |28( as calculated by the JLD(iV) and JLD-RE(A r , N/2) 
procedures using smooth vs non-smooth long-range damping 
functions of Eq. (TJ7J). 



The results presented in the last four rows of Table 
llil show that the present a s -values (first entry for each 
case) agree with the result obtained in Ref. [ll| using the 
asymptotic method to about 10 significant digits. How- 
ever, the JLD-RE(7V, N/2) procedure used here required 
only N = 100 grid points to obtain a s -values for both 
isotopologues with uncertainties of only 0.01%. 

An interesting general point concerns the fact that the 
convergence behavior of the a s -values calculated for Cs2 
and Hc2 at high N was not as smooth as it was for the 
model LJ(2n, n) and realistic Xe-e - potentials; this is 
demonstrated by the contrast between the curves in Fig.[U 
and the 'non-smooth' curves in Fig.[5J We attribute this 
behavior to the discontinuous second and higher deriva- 
tives of the damping function F(x) in the potentials of 
Eqs. (flljf and at the point r = r c . To confirm this 
assertion, the a s calculations for 3 He2 were repeated with 
the Heaviside switching function of Eq. replaced by 
the smooth switching function 




9.8 r Q 



Internuclear distance, r m ; n |au| 



FIG. 9. (Color online) Convergence of calculated scattering 
lengths for the Xe-e", 3,4 He 2 (X 1 E+), and Cs 2 (a 3 Ej) sys- 
tems with respect to the lower bound of the integration inter- 
val, r m in, where ro is the distance where U(r) = 0. 



formation of Eqs. ([50)) and (j3"Tj) formally fixes the lower 
bound for ytg{r) to correspond to r = 0, the exponen- 
tially rapid decay of the wavefunction in the classically 
forbidden region under the short-range repulsive poten- 
tial wall means that in practice that lower bound may be 
set quite a bit closer to the left turning point ro where 
U(r) — 0. To examine this point, Fig. |9] shows the con- 
vergence behavior of calculated values of a s for the three 
real systems of Table |TT] as the lower limit of the integra- 
tion interval is shifted inward from the distance ro. It 
was found that values of r min for which the relative er- 
ror approaches the double-precision numerical-noise limit 
can be defined by the criterion 



H(x) 



1 



-2kx ' 



10 



(49) 



The results obtained in this way, plotted as the 'smooth' 
curves in Fig. [51 clearly confirm our assertion. This 
demonstrates the importance of having potential energy 
functions models which have a very high degree of ana- 
lytic smoothness. 

Our final point concerns the sensitivity of the calcu- 
lated a s -values to the position of the inner boundary 
^min(y = a) of the integration region. While the trans- 



7'0 ) X 



\U(r min )\ > 23 



(50) 



which is based on semi-classical exponential decay of the 
original ^(r) wavefunction in the classically forbidden 
region by a factor of ~ 10~ 10 . For the yt g {r;f,(3) map- 
ping function of Eq. (|31 1) , this implies an inner (left hand) 
boundary point of a = — tan -1 [j3 (r m i n /f — 1)] • 
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IV. CONCLUSIONS 



P 2280. 



This paper presents a very robust and efficient new 
way of calculating scattering lengths. The two-parameter 
mapping function of Eq. ()30|) transforms the conventional 
radial Schrodinger equation ^ into the equivalent form 
of Eq. defined on the finite domain y 6 [a, 1]. For 
arbitrary interaction potentials, neither the solution 4>{y) 
of the transformed equation nor its first derivative <p' '(y) 
are singular anywhere on this interval. As a result, the 
s-wave scattering length a s can be exactly expressed in 
terms of the logarithmic derivative of this transformed 
wave function </>(y = 1) at the right boundary point, as 
specified by Eq. (l33l) . The method does not depend on a 
particular asymptotic form of the potential as long as it 
dies off as fast as 1/r" for n > 4, as r — > oo. 

For well-bound potentials with equilibrium distance r e 
and a limiting (attractive) long-range behavior of l/r n , 
the optimal values of the j/tg( r ) mapping parameters have 
been shown to be f ~ r e and j3 « § — 1, respectively. 
These same mapping parameters also yield efficient so- 
lutions when this approach is used for solving bound- 
state problems. Regardless of the absolute magnitude of 
the scattering length or the number of bound levels sup- 
ported by the potential, the accuracy of the a s calculation 
can easily achieve 10-12 significant digits when working 
in ordinary double-precision arithmetic. It is also shown 
that stable and highly precise values of a s cannot be de- 
termined for analytical potentials which lack high-order 
analytic smoothness throughout the classically allowed 
region. A computer program applying this approach has 
been submitted to the Journal's online data archive [2!| . 

Finally, we note that although it requires a little more 
computational effort to achieve a given level of precision, 
combination of the yoj(r) mapping function of Eq. (pM)) 
with a conventional Numerov wavefunction propagator 
[2l| can yield both accurate a s values and wavefunctions 
4>{y) that can be used for calculating photoassociation 
cross sections and other properties [30] - 
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Appendix: Johnson's log-derivative method 

As was shown in Refs. [23j and [24|, Johnson's quadra- 
ture procedure for outward integration of the Riccati 
equation ([TTj) is based on the two-point finite-difference 
scheme: 



Zfc-l 
1 + Zk-l 



h 2 

— | "'/, "/,■ 



(A.l) 



where £k(Vk) = h 1 z k (y k ), the mesh points are y k 
a + kh, the integration step is h — (b — a)/N, and 



Uk 



Q{Vk) 

Q(y k ) 



k = 0,2,4, 



• ,N 

• ,N- 



(A.2) 



l+(/i 2 /6)Q(y fc 
with weights Wk as in a Simpson quadrature 

1 k = 0,N 
,r, = ^ 4 k = 1,3,5,... ,N-1 

2 k = 2,4,6, ... ,N-2 



(A.3) 



The total number of integration points must be odd, so 
iV must be an even number. 

In the classically forbidden region where Q(yo) < 0, 
the initial value log-derivative solution zq — z(a) can be 
estimated using the semi-classical approximation [l6l[24| : 



z = h 



-Q(yo) - 



Q(yo) 



(A.4) 



It has been verified by numerical calculations [24( 
that the cumulative truncation error of the log-derivative 
method is given by 



= Ch A + 0(h 6 



(A.5) 



As a result, an approximate solution £/t(&) determined 
with fixed integration stepsize (or a fixed number of mesh 
points) hi(Ni) can be extrapolated to zero step size using 
Richardson's formula [25|, [M : 



where 



A 4 - 1 



h2 
hi 



N 2 



(A.6) 



(A.7) 
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